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Abstract 

The requisites for the persistence of small colonies of self-replicating 
molecules living in a two-dimensional lattice are investigated analyt- 
ically in the infinite diffusion or mean-field limit and through Monte 
Carlo simulations in the position- fixed or contact process limit. The 
molecules are modeled by hypercyclic replicators A which are capable 
of replicating via binary fission A + E — > 2 A with production rate s as 
well as via catalytically assisted replication 2A+E — > 3 A with rate c. In 
addition, a molecule can degrade into its source materials E with rate 
7. In the asymptotic regime the system can be characterized by the 
presence (active phase) and the absence (empty phase) of replicators 
in the lattice. In both diffusion regimes, we find that for small values 
of the ratio c/7 these phases are separated by a second-order phase 
transition which is in the universality class of the directed percolation, 
while for small values of s/7 the phase transition is of first order. Fur- 
thermore, we show the suitability of the dynamic Monte Carlo method, 
which is based on the analysis of the spreading behavior of a few active 
cells in the center of an otherwise infinite empty lattice, to address the 
problem of the emergence of replicators. Rather surprisingly, we show 
that this method allows an unambiguous identification of the order of 
the nonequilibrium phase transition. 

PACS: 87.10+e, 87.90.+y, 89.90+n 
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1 Introduction 



The most fundamental event in the history of life was probably the appearance via 
spontaneous creation of a molecule capable of replicating itself (replicator). Given 
a possible mechanism of replication, which in this case is some form of template 
activity, the evolution of such replicators has been extensively investigated through 
the chemical kinetics formalism put forward by Eigen and co-workers in the seven- 
ties (H Those studies have raised a series of objections to the simplistic view 
of the emergence of a complex organism from a collection of competing species of 
replicators. For instance, the finding that the length of a molecule (polynucleotide) 
is limited due to the finite replication accuracy per nucleotide has prompted the 
proposal of models that incorporated some sort of cooperation between the repli- 
cators, such as the molecular catalytic feedback networks termed hypercycles 
These models, however, have attracted their own criticisms since, as pointed out 
by Maynard Smith, giving catalytic support in such molecular networks is in fact 
an altruistic behavior and so they are extremely vulnerable to the presence of par- 
asites, i.e., molecules that do not reciprocate the catalytic support they receive ||. 
A possible solution to this problem is provided by the structured deme formula- 
tion of group selection [Q, where it is assumed that the replicators are spatially 
localized, say, in rock crevices or water droplets, so that the benefits accrued from 
cooperation are directed mostly to the members of the catalytic network ^, 0- 
Yet another successful approach to the problem of resistance against parasites is 
based on a reaction-diffusion system where replication and diffusion taking place on 
an adsorbing surface generate self-organized spiral structures j^, Interestingly, 
as these spatial structures, which greatly increase the stability of the hypercycles 
against parasites, can be viewed as super-organisms that approach is also related 
to the group selection theory Q . 

Since in the prebiotic or chemical evolution context, natural selection is es- 
sentially the dynamics of replicators, it is not surprising that most of the studies 
in this subject have focused almost exclusively on the competition between repli- 
cators, among which the so-called malthusian and hypercyclic replicators are the 
most important ^ . The former corresponds to the simplest reproduction process, 
namely, the binary fission of a parent replicator and is modelled by the chemical 
reaction 

A + E^2A (1) 

where A is the replicator and E is the source materials (mononucleotide resources). 
It is well-known that the concentration of A grows exponentially with the rate con- 
stant s, provided that the concentration of E is kept constant, hence the name 
malthusian replicator. To avoid this explosive growth, one usually imposes a con- 
straint on the total concentration of replicators which can be implemented in prac- 
tice by a dilution flux [Q . Alternatively, one can allow the replicators to be degraded 
by hydrolysis into its mononucleotide components E according to the reaction 

A^E (2) 

which seems a more natural approach to limit the growth of A. 
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As best exemplified by sexual reproduction, there are situations that cannot be 
described by (jl]) since two replicators are necessary to produce a third one. In this 
case the corresponding chemical reaction is 

2A + E A 3A (3) 

which leads to a hyperbolic growth of the concentration of replicator A ||. A 
hypercyclic replicator is defined as one that can replicate itself using both (@) and (|) 
reaction schemes. Actually the term hypercycle derives from the superimposition 
of the catalytic replication cycle (^) on the self-replication cycle (|l|). Of course, the 
limit c = corresponds to the malthusian replicator while s = can be associated 
to an obligatory sexual replicator. 

In contrast to previous works that have concentrated on the competition be- 
tween replicators either of the same kind but with different production rates [jl], |^| or 
of different kinds ||] , in this paper we address a more fundamental problem that has 
received comparatively little attention, namely, the stability of the different kinds 
of replicators, viewed here as an active (ordered) phase of the molecular system, 
against the empty (disordered) phase composed of the resource materials only. This 
lack of interest was probably due to the fact that the usual kinetics formalism used 
to study the dynamics of replicators does not represent the mononucleotide resource 
dynamics explicitly (see however 11 1), thus precluding the study of the issues 



addressed in the present contribution. More pointedly, we consider the dynamics of 
a population of identical hypercyclic replicators on a lattice space both in the de- 
terministic infinite diffusion (mean-field) limit and in the stochastic position-fixed 
(contact process) limit where each replicator on a lattice cell never moves. The 
last limit is particularly interesting because it allows the connection between the 
replicator models and some standard models of noncquilibrium phase transition in 
a lattice (e.g. directed percolation) jlj, |l3|, [uj. As a result, the powerful analytical 
tools of statistical mechanics can be readily used to advance our understanding 
of the evolution of replicators. Of particular relevance is the so-called dynamic 
Monte Carlo method whose idea is to set the system initially in the empty state 
with a seed of replicators in the center of the lattice and then study the subsequent 
spreading of activity [D| [l4| . More importantly, the thorough analysis of both 
limits exposes the limitations of the widely used deterministic chemical kinetics or 
mean-field formalism to study the problem of the emergence of life. 

The remainder of the paper is organized as follows. In Sec. |^ we present the 
set of rules that govern the evolution of a population of hypercyclic replicators in 
a two-dimensional square lattice. The mean-field or infinite diffusion limit, which 
models an ideally mixed medium, is studied analytically in Sec. |[ The results 
are summarized in a phase diagram showing the regions of stability of the empty 
and active regimes in the space of the control parameters of the model. Those 
regions are delimitated by continuous as well as discontinuous transition lines that 
join at a tricritical point. In Sec. |]we study the position-fixed or contact process 
limit using mainly the dynamic Monte Carlo method which allows the computation 
of the critical dynamic exponents that describe quantitatively the spreading of a 
vanishingly small population of replicators. Finally, some concluding remarks are 
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presented in Sec. In particular, we compare the hypercyclic replicator model with 
Schlogl's models of noncquilibrium phase transition in reaction-diffusion systems 

2 The model 

It has been suggested that chemical evolution started with a surface-bonded au- 
tocatalytic chemical network as there are enormous thermodynamics and kinetics 
advantages of surface binding reactions, especially in the case of reactions that re- 
quire unlikely collisions of many reactants fl6}| . The binding must be strong enough 
to keep the reactants on the surface but also flexible enough to allow their slow mi- 
gration on it. Interestingly, if this proposal proves correct it will probably lead to 
the replacement of the popular notion of a primitive soup by that of a primitive 
pizza instead [l7| . Accordingly, we define our replicator model in a two-dimensional 
space consisting of L x L cells in a square toroidal lattice. Each cell is either empty 
or occupied by a replicator and it is assumed that an empty cell contains all source 
material required to assemble a new replicator. The evolution of the population of 
replicators is governed by the following local rules: 

(1) A replicator has a probability 7 of decaying; after decay the cell becomes 
empty. This rule is motivated by the hydrolysis reaction (Q). 

(2) A replicator in one of the four first neighbor cells (von Neumann neighbor- 
hood) of an empty cell can replicate into that cell with probability s. This 
process is referred to as non-catalised self-replication and is motivated by the 
reaction (|l|). 

(3) Regardless of the previous rule, a replicator in the von Neumann neighbor- 
hood of an empty cell can replicate into that cell if there are other replicators 
in the intersection of the Moore neighborhoods of both cells. The probability 
of this type of replication, which is motivated by reaction (^|), is c for each 
pair of replicators. We recall that the Moore neighborhod of a given cell 
consists of its first and second nearest neighbors, adding up to eight cells. 

Hence in the extreme situation where an empty cell is surrounded by eight repli- 
cators, it can become occupied with probability 4s + 16c. To carry out the simu- 
lations we choose the parameters c and s such that 4s + 16c < 1. These rules are 
applied simultaneously to all cells in the lattice so our model can be viewed as a 
two-dimensional stochastic cellular automaton. Actually the model is essentially an 
adaptation to one-membered hypercycles of the spatial cellular automaton model of 
multi-membered hypercycles proposed by Boerlijst and Hogeweg ^|. The dynamics 
defined by the rules given above is manifestly irreversible and, in particular, the 
state characterized by empty cells only is an absorbing state, i.e., a configuration 
from which the system cannot escape. In this sense, the principle of detailed bal- 
ance is broken and the active stationary state is in fact in nonequilibrium. Although 
the more realistic situation is a diffusion-controled reaction where each reactant can 
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move randomly on the lattice, in this paper we choose to study in detail the simpler 
extreme cases of infinite diffusion (mean field) and no diffusion (contact process) . 
Of course, we hope that features common to both limits will be present in the finite 
diffusion situation as well. 



3 The mean-field limit 

The mean-field limit describes exactly an infinite population of reactants in an 
ideally mixed medium and so it is equivalent to the usual chemical kinetics formu- 
lation. Neglecting spatial correlations among cells, i.e., assuming that at any time 
the molecules are distributed randomly over the lattice cells, it is straightforward 
to write the evolution equation for the density of replicators or occupied cells at 
time t, namely, 

pt+i =pt(l-j) + Ap t (1 - p t ) (s + kept) . (4) 

We will consider only the stationary solutions pt+i = pt — p of this equation. The 
absorbing (empty) state p = is always a solution, while the non-zero solutions are 
given by the roots of the quadratic equation 

cp 2 - p(c - s) + 1 - s = (5) 

where we have introduced the dimensionless parameters 

c = and s = — . (o) 

7 7 

This equation has real roots provided that the condition (c + s) 2 > 4c is satisfied. 
In addition we can easily show that: (i) for s < 1 and c > s both roots are positive; 
(ii) for s < 1 and c < s both roots are negative; and (iii) for s > 1 only one of the 
roots is negative. Furthermore, though the solution p = exists in the entire plane 
(c, s), it is stable only if the condition 



dp 



t+i 



Opt 



= l- 7 + 4s<l, (7) 

p=0 



which reduces to s < 1, is satisfied. In the region where the non-zero roots are 
physical (i.e., real and positive), the stable root is always the largest one. In Fig. 
|l] we show the steady-state density of replicators for two different choices of initial 
density. We identify three distinct phases: the absorbing or empty phase (E) 
associated to the solution p — 0; the replicating or active phase (^4) associated to 
the solution p > 0; and the phase labeled (EA) where both solutions p = and 
p > are stable. In this phase, the outcome of the dynamics is not determined by 
the control parameters only but also by the initial abundance of replicators. ^From 
this figure it is clear that the system undergoes a continuous non-equilibrium phase 
transition from phase (E) to phase (^4) at s = s c = 1 and c < 1. Explicitly, near 
this transition the density of replicators can be written as 

(8) 
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Figure 1: Mean- field steady-state concentration of active sites as a function 
of the scaled noncatalized self-replication ratio for (bottom to top) c = 
0, 1, 2, 3 and 4. The initial concentrations are (a) po = 1 and (b) po = .001. 



The continuous transition ends at the tricritical point s t — c t — 1 so that for c > dt 
the transition between phases (E) and (EA), that takes place at s — 2>/d — c, as 
well as the transition between phases (EA) and (A), that occurs at s — 1, are 
discontinuous. In particular, the jumps of the densities of replicators are 

A/?=l-l/g 1/2 c>\ (9) 

at the transition (E) — (EA), and 

Ap = l-l/c 5>1 (10) 

at the transition (EA) — (A). These results are conveniently summarized in the 
phase diagram shown in Fig. ^. It is interesting to note that setting s = c we find, 
close to the tricritical point, 

p^(~s-lf 2 (11) 

so that the tricritical point is not in the same universality class as the transition 
observed in the absence of catalytically assisted replication. 

The interpretation of our results within the prebiotic evolution context leads 
to the conclusion that for finite values of c an obligatory sexual replicator cannot 
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Figure 2: Mean- field phase diagram in the plane (c, s) showing the regions 
of stability of the different steady-state solutions. The continuous transition 
ending at the tricritical point (TCP) is represented as a solid line and the 
broken lines indicate discontinuous phase transitions. 
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emerge spontaneously (i.e., appear at vanishingly small concentrations). For in- 
stance, for s = 0, the minimal initial density of replicators necessary to engender 
a prosperous population is p = 1/2 at the transition point c = 4, vanishing as 
1/c for large c. Actually, an initial colony of replicators is certain to grow from 
vanishingly small concentrations provided that s > 1. These conclusions, however, 
must be taken with caution since in the deterministic limit a vanishingly small con- 
centration means an infinite population of replicators, while one would expect the 
first replicators to show up as a single or a few copies at most. Of course, a proper 
understanding of this emergence phenomenon calls for a stochastic approach, which 
is the subject of the next section. 

4 The position-fixed limit 

The primary aim of this section is to determine what features of the rich phase- 
diagram obtained in the mean-field limit show up also in the opposite limit, where 
the replicators are fixed on the lattice cells. This is a rather challenging enter- 
prise as, at least in the case of equilibrium phase-transitions, there is no totally 
unambiguous way by which one can detect the order of the transition through the 
analysis of finite systems alone [Q. Nevertheless, we tackle this problem using both 
a steady-state approach for finite lattices and the dynamical Monte Carlo method 
for lattices of effectively infinite size. 

First we measure de density p of replicators in the steady state. Our results for 
two system sizes (L — 100 and L = 200) are shown in Figs. || and ^. For each set of 
the control parameters we made runs of 10 5 generations, neglecting the first 2 10 4 
generations and recording p at steps of 200 generations. A generation corresponds 
to the simultaneous update of all lattice cells. Each data point is the arithmetic 
mean of these recorded data. Provided that the population is not extinct, the 
results are independent of the choice of the initial configuration. In particular, we 
have made runs starting from all cells occupied or from a seed of only four clustered 
occupied cells. Furthermore, as in the mean-field limit we have verified that our 
results depend only on the ratios s/j and c/j so throughout the remaining of this 
paper the value of the decay constant is held fixed at 7 = 0.05. These results 
clearly indicate the existence of a phase transition separating the active (p > 0) 
and empty (p = 0) phases of the model. Assuming that in the neighborhood of 
the transition points the density of replicators goes to zero as p ~ (s — Sc) 13 and 
using the least-square method we can estimate both the critical replication rate s c 
and the critical exponent (5 > 0. These estimates are presented in table [j] and the 
quality of the fitting can be appreciated from Fig. ||. The statistical errors are of 
order of 10 -3 but the systematic errors, which are due mainly to the difficulty to 
carry out long runs very close to the transition point, are probably much larger. 
We note that our estimates of (3 for small values of c indicate that in this regime the 
replicator model belongs to the so-called (2 + 1) directed-percolation universality 
class for which j3 = 0.59 ± 0.02 It should be pointed out that in finite systems 
the active regime is a meta-stable one as there is always a finite probability that 
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Figure 3: Average density of replicators p as a function of s for (left to right) 
c = 5.70,5.13,3.94,2.58,0.99 and 0. The lattice sizes are L = 100(x) and 
200(0). 



the colony becomes extinct due to fluctuations in the stochastic dynamics. Since 
this extinction probability increases towards 1 as the critical point is approached, 
it is very difficult to obtain reliable estimates of s c and (3 by means of numerical 
simulations in the steady-state regime. Of course we are aware that if the transition 
happens to become discontinuous at some value of c > then the assumption of a 
power-law singularity at the transition point breaks down (actually /3 = in this 
case). In fact, the (anomalous) continuous decrease of (3 as c increases (see table 
is an indication that this might be the case. Our careless use of the power-law 
assumption is intentional and aims to illustrate the difficulty of detecting the order 
of a phase transition using results of steady-state simulations: the abrupt variation 
of p at criticality observed in Fig. || for certain values of c, which might indicate the 
occurrence of a first-order irreversible transition, can be explained as a continuous 
transition with a small exponent f3 as well. 

We turn now to the analysis of the spreading behavior of a small colony of 
replicators settled initially in the center of an otherwise empty lattice of infinite 
size. More pointedly, the initial colony is composed of four replicators located in the 
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Figure 4: Logarithm plot of the average density p as a function of (s — s c ) for 
(top to bottom) c = 5.70, 5.13, 3.94, 2.58, 0.99 and 0. The values of s c = s c (c) 
for the different choices of c are given in table [l] and the straight lines are 
the numerical fitting obtained with those data. Only the data for L = 200 
are presented. 



c s c j3 

1.625 061 

0.992 1.400 0.60 

2.576 1.004 0.47 

3.936 0.608 0.30 

5.128 0.210 0.18 

5.7008 0.004 0.12 



Table 1: Estimates of the critical point s c and the critical exponent j3 as- 
suming the power-law singularity p ~ (s — s c )^. 
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von Neumann neighborhod of the central empty cell. Finite size effects are absent 
because the lattice size is taken large enough so that during the time we follow the 
evolution of the colony the replicators can never reach the lattice boundaries. This 
of course sets an upper limit to the number of generations we can follow the colony 
and so, in particular, we let the population evolve up to typically t — 10 4 . As usual, 
we concentrate on the time dependence of the following key quantities |Q : (i) the 
average number of replicators n(t); (ii) the survival probability of the colony P(t); 
and (iii) the average mean-square distance over which the replicators have spread 
R 2 (t). For each time t we carry out M = 2 10 5 independent runs, all starting with 
the same initial colony. Hence P(t) is simply the fraction of runs for which there is 
at least one replicator in the lattice at time t. Since n(t) is an average taken over all 
runs including those which have already been extinct at generation t, the average 
number of replicators per surviving run is given by the ratio N(t) = n(t)/P(t). 
Furthermore, noting that R 2 (t) is averaged only over the surviving runs, we can 
define the fractal dimension df of the surviving colonies of replicators at a given 
time t as N ~ R d f . 

At the transition points we expect that the measured quantities obey the fol- 
lowing scaling laws Jl2) 

n(t) ~ f> (12) 
P(t) ~ r s (13) 
R 2 (t) ~ t z (14) 

where S, rj and z are dynamic exponents which are related with the fractal dimension 
of the clusters of replicators through the equation 

d f = 2*±*. (15) 

z 

In principle, these scaling laws are valid for continuous as well as discontinuous 
phase transitions, though the scaling relations between the exponents, such as the 
'hyperscaling' relation jlS 

]^dz-ri = 25 (16) 

where d is the lattice dimension, hold only in the case of continuous transitions. 

In Figs. [| and || we present log-log plots of n{t) and P(t), respectively, as 
functions of t in the vicinity of the critical point for 5=0. The dependence of 
R 2 {t) on t is not shown since, near criticality, the curves for different values of s 
are clustered together and do not reveal any qualitatively relevant information on 
the colony evolution. The asymptotic straight lines observed in these figures are 
the signature of critical behavior while upward and downward deviations indicate 
supercritical and subcritical behaviors, respectively. A precise estimate for the 
critical exponents is obtained by considering the local slopes of the curves shown 
in the previous figures. For instance, the local slope 6(t) is defined by (l4[ ^U) 

_ m = m[P(t)/P(t/B)] 

m 5 
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Figure 5: The log-log plot of n(t) as a function of t for c = and (top to 
bottom) a = 1.80, 1.632, 1.630, 1.629, 1.628, 1.627, 1.626, and 1.52. 



which for large t behaves as 

5(t)~5 + - t (18) 

where a is a constant. Analogous expressions hold for rj(t) and z(t). Hence plots of 
the local slopes as functions of 1/t allow the calculation of the critical exponents. 
Using this procedure we find s c — 1.628 ± 0.001 which yields the exponents r\ = 
0.23 ± 0.01, S = 0.45 ± 0.01 and z = 1.13 ± 0.01. The error in S c is estimated 
by determining two values of s as close as possible to the critical point for which 
upward and downward deviations can be observed, while the errors in the critical 
exponents are, as usual, the statistical errors obtained by fitting the local slopes 
by straigh lines in the large t regime. Our exponents are in good agreement with 
those of the (2 + 1) directed-percolation |Q and satisfy very well the hyperscaling 
relation ( |l6| ) thus indicating that the transition in the limit c = is continuous, as 
expected. 

The results of the spreading analysis for the other extreme case, 5 = 0, which 
models a population of obligatory sexual replicators are shown in Figs. ^ and |[ 
Although the dependence of lnP(i) on t is similar to that observed in the previ- 
ous case, the behavior pattern of lnn(£) (see Fig. ^) is rather different: in the 
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Figure 6: Same as fig. || but for P(t). 
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supercritical regime (i.e., c > c c ) n(t) first increases reaching a maximum, then 
decreases reaching a minimum and finally starts to increase monotonically again. 
This pattern is illustrated better in Fig. ||, which shows the spreading results for 
s = 0.2. Analysis of these figures, which exceptionally show the colony evolution 
up to 2 10 4 generations, suggests that a flat line separates the supercritical and 
the subcritical regimes implying thus the vanishing of the exponent 77 at the tran- 
sition point. Furthermore, the qualitatively distinct behavior patterns of hm(t) 
illustrated in Figs. |E|, [?] and ^ can be used to identify unambiguously the order 
of the nonequilibrium phase transition and hence to estimate the location of the 
tricritical point. To appreciate how the time dependence of lnn(t) in the supercrit- 
ical regime changes continuously from the simple monotonic increase for s = to 
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the complicated behavior described above for c = we present in Figs. 10 and 
log-log plots of n(t) as a function of t for s = 0.6 and s = 1.0, respectively. In fact, 
analysis of Figs. to [ll] suggests that the turning point between those distinct be- 
havior patterns occurs when the maximum and the minimum of n(t) coincide, i.e., 
the critical curve lnn(i) vs. Int has an inflection point (see Fig. |l(j). The value of 
s and the corresponding c c at which this behavior occurs are then identified as the 
coordinates of the tricritical point. Applying this procedure we find St — 0.60±0.04 
and c t = 4.00 ± 0.20. 

The transition points and the dynamic exponents obtained via the scaling laws 
([12])- ([ll]) and via the analysis of the local slopes are summarized in table ^. The 
errors in the estimates of the transition points are calculated as described before 
for the case c = 0. Except for that case, we refrain from presenting the (statistical) 
errors in the exponents since the systematic errors are unusually large, due probably 
to the crossover behavior among (at least) three different universality classes. For 
instance, in the vicinity of the transition point for s = 0.6, analysis of the local slope 
S(t) up to t — 500 indicates a clear tendency to the asymptotic value 6 « 1 while 
for t > 500 the tendency suddenly changes towards the asymptotic value S » 0.6. A 
similar phenomenon occurs for s = 1.0 also: the initial tendency is towards S w 0.6 
and then changes towards 8 « 0.45 for larger times. As a result the estimate of 
the exponents becomes strongly dependent on the precise location of the transition 
points, which requires even better statistics as well as much longer runs. We leave 
this interesting research line which includes, for instance, the identification of the 
universality class of the tricritical point to a future, more technical contribution. 

The evidences in support of our claim that, similarly to the mean-field limit, 
in the position-fixed limit the nonequilibrium phase transition between the empty 
and active phases is discontinuous for small s and so there is a tricritical point in 
the phase diagram of the model are threefold: First, the vanishing of the exponent 
77 in this range signalizes a distinct asymptotic behavior of the average number of 
replicators n(t). Second, the hyperscaling relation ( |l6| ) is clearly violated for small 
s while it is satisfied in the regime where we expect the transition to be a second- 
order one. Third, using Eq. ( |l5| ) and the data of table |^ we find df w d = 2 in the 
region of small s indicating that the clusters of replicators are not fractal objects, 
in contrast to the clusters observed in the vicinity of a second-order transition (we 



find, for instance, df « 1.21 for c = 0). This point is illustrated in Figs. O and 13 
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Figure 7: The log-log plot of n(t) as a function of t for s = and (top to 
bottom) c = 5.760, 5.712, 5.709, 5.704, 5.702, 5.701, 5.680 and 5.664. 
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Table 2: Critical dynamic exponents calculated from the slopes of the 
straight lines at the transition points c c (s). 
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Figure 8: Same as fig. |^ but for P(t). 
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Figure 9: The log-log plot of n(t) as a function of t for s = 0.2 and (top to 
bottom) c = 5.20, 5.168, 5.160, 5.152, 5.144, 5.136 and 5.120. 
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Figure 11: The log-log plot of n(t) as a function of t for s = 1.0 and (top to 
bottom) c = 2.720, 2.600, 2.592, 2.584, 2.576, 2.568, 2.552 and 2.544. 
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Figure 12: Snapshot of the lattice configuration at t = 10 4 showing the 
colony of replicators (dots). The parameters are c = and s = 1.69. The 
initial colony of four replicators was placed in the middle of the 200 x 200 
lattice. 



which show snapshots of typical colonies at t = 10 for the two extreme cases. In 
both figures the relative distances to the critical points are the same. The reason for 
the colonies of obligatory replicators to be much denser than those of malthusian 
replicators is that the number of replicators per surviving colony (N(t) ~ t n+s ) 
increases roughly as t for the former and as t - 7 for the latter. In addition, the 
average square distance over which the initial colony has spread from the center of 
the lattice at time t is roughly of the same order in both extremes as indicated by 
the values of the exponent z. 

Finally, in Fig. |l4| we present the phase diagram for the position-fixed limit. 
In this case there are only two phases, namely, the empty phase (E) characterized 
by a vanishing probability of survival Poo = lim^oo P(t) — and the phase (EA) 
where the active and empty states can occur with probabilities Poo and 1 — Poc 
(see Figs. || and ||), respectively. We note that the population size is effectively 
unlimited so that extinction is not certain to occur as in the steady-state analysis 
of finite systems. It is worth emphasizing that for finite production rates one has 
Poo < 1 and so there is always a nonvanishing probability of extinction. (Actually, 
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Figure 13: Same as Fig. 12 but for s = and c = 5.92. 
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Figure 14: Phase diagram in the plane (c, s) for the position-fixed limit. 
The continuous and discontinuous transition points are represented by the 
symbols O an d A, respectively. The error bars represent the uncertainty 
in the location of the tricritical point (TCP) while for the other data points 
the error bars are smaller than the symbol sizes. For sake of comparison the 
mean-field phase-diagram (solid and broken lines) is also presented. 



Poo depends on the number and location of the replicators in the initial colony, 
but the conditions we have chosen are the most relevant for the emergence of life 
problem.) This result contrasts to our findings in the mean-field limit that the 
fixed point associated to the empty phase becomes unstable for s > 1 so that 
even starting with a vanishingly small concentration of replicators the population 
never dies out. As mentioned before, the reason for this discrepancy is not the 
difference in mobility of the replicators but the fact that the mean-field analysis 
actually considers an infinite population and hence it fails to take into account the 
stochastic fluctuations that could drive a small population to extinction. 

For sake of completeness, we should mention that we have also carried out 
a similar analysis for one-dimensional lattices (chains). While the results for the 
mean-field limit are of course the same (provided we properly redefine the dimcn- 
sionless parameters c and s), the fixed-point limit has some distinct features that 
are worth-mentioning. In particular, we find no evidence for a first-order tran- 
sition; instead, we find that in both extremes c = and s = the empty and 
active phases are separated by a second-order phase transition that belongs to the 
(1 + 1) directed-percolation universality class fl2| . Interestingly, in one-dimension 



22 



the steady-state analysis of finite chains is rendered practically useless by the very 
pronounced finite-size effects, which are probably due to the proximity to the lower 
critical dimension of the model. 

5 Conclusion 

In this paper we have focussed on the prior step in the evolution of life: What are 
the necessary conditions for small colonies of molecules capable of making copies 
of themselves via some template mechanism to persist? This step must be passed 
before one can consider issues such as the outcome of the competition between the 
replicators and their defective copies ||, or between different kinds of replicators 
. As a model of replicator we have considered the well-established hypercyclic 
replicator (one-membered hypercycle) which incorporates two independent mech- 
anisms of replication, namely, the direct template replication reaction ([!]) and the 
catalytically-assisted template replication reaction (||), whose rates are proportional 
to the parameters s and c, respectively. Furthermore, motivated by the modern the- 
ories for the evolution of life that suggest a scenario of diffusion-controlled chemical 
reactions taking place on adsorbing surfaces (probably pyrite) where each reactant 
can move randomly on the surface |l6|, [lTj], we have considered a two-dimensional 
lattice model where each reactant can occupy one of the lattice cells. 

Since the diffusion process of reactants complicates considerably the analysis, 
we have focussed on the two extreme situations: the infinite diffusion or mean- 
field limit and the position- fixed or contact process limit. The expectation is that 
features common to both limits should also be present in the more realistic, inter- 
mediate situations. In both cases we found rich phase diagrams showing the regions 
in the plane (c, s) where the replicators persist (active phase) and die out (empty 
phase): these regions are separated by second-order nonequilibrium phase transi- 
tions which turn into first-order ones at tricritical points. The dynamic Monte Carlo 
method has proven very well suited to our investigation of the position-fixed limit 
not only because the method is based on the analysis of the spreading behavior of 
a small colony of active cells, which is exactly the problem we are interested in, but 
because, rather surprisingly, it allows an unambiguous identification of the order 
of the nonequilibrium phase transition. In addition, we show that the continuous 
transition is in the universality class of the (2 + I) directed percolation. 

Some remarks on the apparent similarity between Schlogl's first and second 
models and the replicator models studied in this paper are in order Jlq] . In fact, 
irreversible versions of Schlogl's first and second models are recovered when the 
reactant E is eliminated from reactions (Q) and (||), respectively, so that the exis- 
tence of an empty cell containing source materials is not required for replication. 
Interestingly, this difference is not important in the case of malthusian replicators 
(c = 0) since this model has the same critical behavior as Schlogl's first model, 
namely, a second-order phase transition which is in the same universality class of 
the (d + 1) directed percolation (the mean- field limit is obtained for d > 4) pl[ . 
The comparison between Schlogl's second model and obligatory sexual replicators 
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(s = 0) is more involved. On the one hand, the mean-field analysis of Schlogl's sec- 
ond model predicts a first-order phase transition |f5| but Monte Carlo calculations 
indicate that the transition for d < 4 is a second-order one which is actually in the 
same universality class as the transition in Schlogl's first model . On the other 
hand, our results for the case s — show that a directed-percolation-like, second- 
order phase transition takes place for d = 1 only, the transition being discontinuous 
for higher dimensions (d > 2). We should mention, however, that the Monte Carlo 
implementation of Schlogl's second model actually allows the diffusion of reactants 
to neighboring cells and, in addition, allows a cell to shelter more than one reactant 
[ p2[ so that a comparison with the position-fixed limit may not be appropriate. 
In any event, it is our opinion that the hypercyclic replicator model should not be 
viewed as a mere variant of Schlogl's models; rather it is a well-established model of 
chemical evolution 0, [|, [j| which, as far as we know, has not been studied beyond 
the mean-field limit. 

Although the simple replicator model considered in this paper turned out to 
be a quite exciting model of noncquilibrium phase transitions, we should not lose 
sight of the original purpose of this work and so, at this stage, it is important to 
highlight the relevance of our results to the origin of life issue. In fact, the mere 
existence of a phase transition between the empty and the active regimes poses a 
difficulty to our scenario of the emergence of life since the scaled production rate 
s of the spontaneously created self-replicating molecule must be larger than some 
threshold value already at the outset. Though increasing the mobility of the reac- 
tants decreases this threshold somewhat, that scenario would be more plausible if 
replicators with vanishingly small production rates could also thrive. The situation 
becomes even worse in the case of first-order transitions: in the deterministic mean- 
field limit the initial abundance of the spontaneously created replicators should be 
large as well, while in the stochastic position-fixed limit the probability of survival 
in the vicinity of the transition point is some orders of magnitude smaller than in 
the case of a second-order transition (see Figs. |6|and||). Furthermore, our results 
indicate that some important conclusions, such as the certainty of survival for s > 1 
or the role played by the initial concentration of replicators near a discontinuous 
transition, are actually artifacts of the deterministic formalism commonly used to 
study chemical evolution. 

In summary, our results show the necessity of adding some new elements to 
the standard scenario for the emergence of life whose effect would be to avoid the 
phase-transitions, allowing thus inefficient replicators to thrive at this first stage 
of life. Only then one can invoke natural selection and imperfect replication to 
boost the replication rates. In addition, our results point to the inadequacy of the 
deterministic mean-field or chemical kinetics formulation to address the origin of 
life issue and suggest as an alternative stochastic formulation the dynamic Monte 
Carlo method that has been extensively used in the physics literature to study 
nonequilibrium phase transitions. 
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